getwd()
setwd(getwd())
install.packages(c("foreign", "stringr", "gdata", "car","zoo","tidyr","RColorBrewer","dplyr", "ggplot2", "readxl", "readr","WDI", "countrycode"))
BasePackages <- c("foreign", "stringr", "gdata", "car", "zoo", "tidyr", "RColorBrewer", "dplyr", "ggplot2", "readxl", "readr","gdxrrw")
#setwd("D:\thasegawa\ExtremeEvent\R")
library(plyr)
library(reshape2)
library(ggplot2)
library(plyr)
library(grid)
library(scales)
library(gdxrrw)
library(dplyr)
colnames <- c("VAR","Year","Region","SSP","GCM","RCP","CO2fert","SYR","N","Crop","Value")
colnamesNocc <- c("VAR","Year","Region","SSP","Crop","Value")

regionlabel <- function(y){
  levels(y$Region)[levels(y$Region)=="WLD"] <- "World"
  levels(y$Region)[levels(y$Region)=="USA"] <- "United States"
  levels(y$Region)[levels(y$Region)=="XE25"] <- "EU25"
  levels(y$Region)[levels(y$Region)=="XER"] <- "Rest of Europe"
  levels(y$Region)[levels(y$Region)=="TUR"] <- "Turkey"
  levels(y$Region)[levels(y$Region)=="XOC"] <- "Oceania"
  levels(y$Region)[levels(y$Region)=="CHN"] <- "China"
  levels(y$Region)[levels(y$Region)=="IND"] <- "India"
  levels(y$Region)[levels(y$Region)=="JPN"] <- "Japan"
  levels(y$Region)[levels(y$Region)=="XSE"] <- "Southeast Asia"
  levels(y$Region)[levels(y$Region)=="XSA"] <- "Rest of South Asia"
  levels(y$Region)[levels(y$Region)=="CAN"] <- "Canada"
  levels(y$Region)[levels(y$Region)=="BRA"] <- "Brazil"
  levels(y$Region)[levels(y$Region)=="XLM"] <- "Rest of South America"
  levels(y$Region)[levels(y$Region)=="CIS"] <- "FSU"
  levels(y$Region)[levels(y$Region)=="XME"] <- "Middle East"
  levels(y$Region)[levels(y$Region)=="XNF"] <- "North Africa"
  levels(y$Region)[levels(y$Region)=="XAF"] <- "Sub-Sa Africa"
  z <- y
  return(z)
}

croplabel <- function(y){
  levels(y$Crop)[levels(y$Crop)=="PDR"] <- "rice"
  levels(y$Crop)[levels(y$Crop)=="WHT"] <- "wheat"
  levels(y$Crop)[levels(y$Crop)=="GRO"] <- "other grains"
  levels(y$Crop)[levels(y$Crop)=="OSD"] <- "oil crops"
  z <- y
  return(z)
}
fsize <- 20
MyThemeLine <- theme_grey() +
  theme(
    panel.border=element_rect(fill=NA),
    title=element_text(size=fsize,face="bold"),
    axis.title=element_text(size=fsize,face="bold"),
    axis.text = element_text(hjust=1,size = fsize, angle = 0),
    axis.line=element_line(colour="black"),
    panel.background=element_rect(fill = "white"),
    #    panel.background=element_blank(),
    #    panel.grid.major=element_line(linetype="dashed",colour="grey",size=0.5),
    panel.grid.major=element_blank(),
    strip.background=element_rect(fill="white", colour="white"),
    strip.text.x = element_text(size=fsize, colour = "black", angle = 0,face="bold"),
    legend.title = element_text(size=fsize),
    legend.text = element_text(size=fsize)
  )

#Defining function
#Making two figures
twofigure <- function(p1,p2,x){
  grid.newpage()
  l <- grid.layout( nrow=1, ncol=2, heights=unit(c(1, 1.25),c("null","null")), widths=unit(c(1, 1),c("null","null")))
  png(filename=x, width=2400, height=1500,res=200)
  v <- viewport(layout = l) # [l]を作図領域として[v]に格納する
  pushViewport(v) #グラフを埋め込む前に，pushViewport()関数に作図領域を指定する。
  print(p1, vp = viewport(layout.pos.row=1,layout.pos.col=1))
  print(p2, vp = viewport(layout.pos.row=1,layout.pos.col=2))
  popViewport()  #作図を終了する
  dev.off()  
  return()
}

fourfigure <- function(p1,p2,p3,p4,x){
  grid.newpage()
  l <- grid.layout( nrow=2, ncol=2, heights=unit(c(1, 1.25),c("null","null")), widths=unit(c(1, 1),c("null","null")))
  png(filename=x, width=2400, height=1500,res=200)
  v <- viewport(layout = l) # [l]を作図領域として[v]に格納する
  pushViewport(v) #グラフを埋め込む前に，pushViewport()関数に作図領域を指定する。
  print(p1, vp = viewport(layout.pos.row=1,layout.pos.col=1))
  print(p2, vp = viewport(layout.pos.row=2,layout.pos.col=1))
  print(p3, vp = viewport(layout.pos.row=1,layout.pos.col=2))
  print(p4, vp = viewport(layout.pos.row=2,layout.pos.col=2))
  popViewport()  #作図を終了する
  dev.off()  
  return()
}
eightfigure <- function(p1,p2,p3,p4,p5,p6,p7,p8,x){
  grid.newpage()
  l <- grid.layout( nrow=2, ncol=4, heights=unit(c(1, 1.25),c("null","null")), widths=unit(c(1, 1),c("null","null")))
  png(filename=x, width=2400, height=1500,res=200)
  v <- viewport(layout = l) # [l]を作図領域として[v]に格納する
  pushViewport(v) #グラフを埋め込む前に，pushViewport()関数に作図領域を指定する。
  print(p1, vp = viewport(layout.pos.row=1,layout.pos.col=1))
  print(p2, vp = viewport(layout.pos.row=2,layout.pos.col=1))
  print(p3, vp = viewport(layout.pos.row=1,layout.pos.col=2))
  print(p4, vp = viewport(layout.pos.row=2,layout.pos.col=2))
  print(p5, vp = viewport(layout.pos.row=1,layout.pos.col=3))
  print(p6, vp = viewport(layout.pos.row=2,layout.pos.col=3))
  print(p7, vp = viewport(layout.pos.row=1,layout.pos.col=4))
  print(p8, vp = viewport(layout.pos.row=2,layout.pos.col=4))
  popViewport()  #作図を終了する
  dev.off()  
  return()
}
#colnames <- c("VAR","Year","Region","SSP","GCM","RCP","CO2fert","SYR","N","Crop","Value")
Dataextract <- function(v,reg,cr,yr){
  xhist <- pdata[pdata$VAR==v & pdata$Year==yr & pdata$Region==reg & pdata$SSP=="SSP2" & pdata$Crop==cr, c(-1,-2,-3,-4,-10)]
  names(xhist) <- c("GCM","RCP","CO2fert","SYR","N","Value")	#Rename the column
  nocc <- NoCC_load[NoCC_load$VAR==v & NoCC_load$Year==yr & NoCC_load$Region==reg & NoCC_load$SSP=="SSP2" & NoCC_load$Crop==cr, c(6)]  
  minv <- min(min(xhist[6]),nocc)
  if (v=="PPOP_UNSH"){
    maxv <- max(xhist[6],nocc)
    bw <- maxv/100
  } else {
    maxv <- max(xhist[6],nocc)
    bw <- diff(range(xhist$Value))/100
  }
  mainlabel <- paste(reg,cr,yr)
  interact<-interaction(xhist$RCP, xhist$CO2fert, sep=" : ")
  mu <- ddply(xhist, "interact", summarise, grp.mean=mean(Value))
  #  v5 <- ddply(xhist,"interact",summarise, grp.quantile=quantile(Value,0.05))
  px <- ggplot(data=xhist, aes(x=Value,fill=interact)) +
    #    geom_density(alpha = 0.3, adjust = 0.2, position = "identity", stat = "density") +
    geom_histogram(alpha = 0.3, position = "identity", aes(xhist=..density..), binwidth=bw, color="black") +
    xlim(minv-abs(minv)*0.01, maxv+abs(maxv)*0.01) +
    ggtitle(mainlabel) +
    ylab("density") + xlab(labelname) +
    MyThemeLine +
    geom_vline(data=mu,aes(xintercept=grp.mean, color=interact), linetype="dashed", size=1)
  #  geom_vline(data=v5,aes(xintercept=grp.quantile, color=interact), linetype="solid", size=1) +
  #  geom_text(data=mu,aes(x=grp.mean,y=0.023,label=grp.mean),angle=-90,size=5,vjust=-bw/10,hjust=-0.0001) +
  # geom_text(data=v5,aes(x=grp.quantile,y=0.023,label=grp.quantile),angle=-90,size=5,vjust=-bw/10,hjust=-0.0001) +
  #  geom_vline(aes(xintercept=nocc), color="black", size=1.2) 
  #  geom_text(aes(x=nocc1,y=0.023,label=nocc1),angle=-90,size=5,vjust=-bw/10,hjust=-0.0001)
  #  geom_area(mapping = aex(y = ifelse(y<0.2)))
  
}


#---- making directory
dirname0 <- paste('../output/png/', sep = '')
dir.create(dirname0)

var.in <- c("PCAL","PPOP_UNSH","YEXO","XPRP")

for(k in 1:length(var.in)){
  varname <- var.in[k]
  
  #  dirname <- paste('../output/png/', varname, sep = '')
  dirname <- paste(dirname0, varname, sep = '')
  dir.create(dirname)
  
}
#--- data load
pdata <- rgdx.param('../data/Dataall_202050.gdx','Data')
names(pdata) <- colnames
NoCC_load <- rgdx.param('../data/Dataall_202050.gdx','Data_Nocc')
names(NoCC_load) <- colnamesNocc
yearname <- c('2050')

pdata$RCP <- factor(pdata$RCP, levels=c("RCP85","RCP26"))
pdata$Region <- factor(pdata$Region, levels=c("USA","CAN","JPN","TUR","XER","XE25","XOC","CIS","XNF","CHN","BRA","IND","XLM","XSE","XME","XSA","XAF","WLD"))

fao <- rgdx.param('../data/Risk_fao.gdx','risk_fao_tot')
fao <- rgdx.param('../data/Risk_fao.gdx','pcal_fao_tot')
names(fao) <- c("Reg","Year","Value")
faoworld <- fao %>%
  filter(Reg=="WLD")
#            

#-----PPOP_UNSH & PCAL: Histogram for a single region ----
var.in <- c("PCAL")
var.in <- c("PPOP_UNSH")
var.in <- c("PCAL","PPOP_UNSH","XPRP","YEXO")
regionlist <- c("WLD","IND","XSA","XAF","BRA")
regionlist <- c("WLD")
crop_list <- c("PDR","WHT","GRO","OSD") 
nocc <- 0
yearlist <- c("2020","2050")
yearlist <- c("2050")
cropname <- c("dum")



for(k in 1:length(var.in)){
  varname <- var.in[k]
  if (varname=="PCAL"){
    labelname <- c("calorie intake [kcal/person/day]")
    cropnum <- 1
  } else if (varname=="PPOP_UNSH"){
    labelname <- c("population at risk of hunger [billion]")
    cropnum <- 1
  } else if (varname=="YEXO"){
    labelname <- c("yield [tonne/ha]")
    cropnum <- 4
  } else if (varname=="XPRP"){
    labelname <- c("Price index [-]")
    cropnum <- 4
  } else {
    labelname <- c("aaaaaaaa")
  }
  for(r in 1:length(regionlist)){
#    yearname <- yearlist[y]  
    regionname <- regionlist[r]
    
    
    if (cropnum==4){
      for (y in 1:length(yearlist)){
        filenamehist <- paste('../output/png/', varname ,'/',regionname,'_',yearname,'.png', sep = '')
        cropname <- crop_list[1]
        nocc <- NoCC_load[NoCC_load$VAR==varname & NoCC_load$Year==yearname & NoCC_load$Region==regionname & NoCC_load$SSP=="SSP2" & NoCC_load$Crop==cropname, c(6)]  
        plot1 <- Dataextract(varname,regionname,cropname,yearname)
        cropname <- crop_list[2]
        nocc <- NoCC_load[NoCC_load$VAR==varname & NoCC_load$Year==yearname & NoCC_load$Region==regionname & NoCC_load$SSP=="SSP2" & NoCC_load$Crop==cropname, c(6)]  
        plot2 <- Dataextract(varname,regionname,cropname,yearname)
        cropname <- crop_list[3]
        nocc <- NoCC_load[NoCC_load$VAR==varname & NoCC_load$Year==yearname & NoCC_load$Region==regionname & NoCC_load$SSP=="SSP2" & NoCC_load$Crop==cropname, c(6)]  
        plot3 <- Dataextract(varname,regionname,cropname,yearname)
        cropname <- crop_list[4]
        nocc <- NoCC_load[NoCC_load$VAR==varname & NoCC_load$Year==yearname & NoCC_load$Region==regionname & NoCC_load$SSP=="SSP2" & NoCC_load$Crop==cropname, c(6)]  
        plot4 <- Dataextract(varname,regionname,cropname,yearname)
        fourfigure(plot1,plot2,plot3,plot4,filenamehist)
      }
    }else if (cropnum==1){   
      filenamehist <- paste('../output/png/', varname ,'/',regionname,'.png', sep = '')
      cropname <- c("dum")
      yearname <- yearlist[1]
      nocc <- NoCC_load[NoCC_load$VAR==varname & NoCC_load$Year==yearname & NoCC_load$Region==regionname & NoCC_load$SSP=="SSP2" & NoCC_load$Crop==cropname, c(6)]  
      Dataextract(varname,regionname,cropname,yearname)
      plot1 <- px
      yearname <- yearlist[2]
      nocc <- NoCC_load[NoCC_load$VAR==varname & NoCC_load$Year==yearname & NoCC_load$Region==regionname & NoCC_load$SSP=="SSP2" & NoCC_load$Crop==cropname, c(6)]  
      plot2 <- Dataextract(varname,regionname,cropname,yearname)
      twofigure(plot1,plot2,filenamehist)
      ggsave(plot1,filename=filenamehist, dpi=250,width=8, height=6)
      
      
    }
  }
}
#----End --- PPOP_UNSH & PCAL: Histogram for a single region----

# Figure1
#-----PPOP_UNSH & PCAL: Ribbon for a single region ----
var.in <- c("PCAL")
var.in <- c("PPOP_UNSH")
regionlist <- c("WLD","IND","XSA","XAF","BRA")
regionlist <- c("WLD")

crop_list <- c("PDR","WHT","GRO","OSD") 
nocc <- 0
cropname <- c("dum")
yearlist <- c("2020","2050")
ssp_list <- c("SSP2","SSP3")
x <- pdata
y <- NoCC_load

for(k in 1:length(var.in)){
  varname <- var.in[k]
  if (varname=="PCAL"){
    labelname <- c("calorie intake [kcal/person/day]")
    cropnum <- 1
  } else if (varname=="PPOP_UNSH"){
    labelname <- c("population at risk of hunger [billion]")
    cropnum <- 1
  } else if (varname=="YEXO"){
    labelname <- c("yield [tonne/ha]")
    cropnum <- 4
  } else if (varname=="XPRP"){
    labelname <- c("Price index [-]")
    cropnum <- 4
  } else {
    labelname <- c("aaaaaaaa")
  }
  for(r in 1:length(regionlist)){
    regionname <- regionlist[r]
    
    if (cropnum==4){
      crop_list <- c("PDR","WHT","GRO","OSD") 
    }else if (cropnum==1){   
      crop_list <- c("dum")
    }
    for (crn in 1:length(crop_list)){
      cropname <- crop_list[crn]
      #        for (ssp in 1:length(ssp_list)){
      flag <- 0
      xrib <- x[x$VAR==varname & x$RCP=="RCP85" & x$Region==regionname & x$Crop==cropname, c(-1,-3,-6,-10)]
      #            nocc1 <- y[y$VAR==varname & y$Region==regionname & y$SSP=="SSP2" & y$Crop==cropname, c(6)]  
      names(xrib) <- c("Year","SSP","GCM","XXX","SYR","N","Value")	#Rename the column
      
      minv <- min(min(xrib[7]))
      if (varname=="PPOP_UNSH"){
        maxv <- max(max(xrib[7]))
        bw <- maxv/100
      } else {
        maxv <- max(max(xrib[7]))
        bw <- diff(range(xrib$Value))/100
      }
      
      mainlabel <- paste(regionname,cropname)
      faoworld2 <- faoworld
      faoworld2[2] <- lapply(faoworld2[2],as.integer)
      faoworld2[2] <- faoworld2[2]+1990
#      interact<-interaction(xrib$RCP, xrib$SSP, sep=" : ")
      
      xrib1 <- ddply(xrib,c("Year","SSP"),transform,max=max(Value))
      xrib2 <- ddply(xrib1,c("Year","SSP"),transform,min=min(Value))
      xrib5 <- ddply(xrib2,c("Year","SSP"),transform,v1=quantile(Value,0.01))
      xrib6 <- ddply(xrib5,c("Year","SSP"),transform,v2=quantile(Value,0.02))
      xrib7 <- ddply(xrib6,c("Year","SSP"),transform,v5=quantile(Value,0.05))
      xrib8 <- ddply(xrib7,c("Year","SSP"),transform,v10=quantile(Value,0.1))
      xrib9 <- ddply(xrib8,c("Year","SSP"),transform,v20=quantile(Value,0.2))
      xrib10 <- ddply(xrib9,c("Year","SSP"),transform,v33=quantile(Value,0.3333))
      xrib105 <- ddply(xrib10,c("Year","SSP"),transform,v50=quantile(Value,0.5))
      xrib11 <- ddply(xrib105,c("Year","SSP"),transform,v66=quantile(Value,0.6666))
      xrib12 <- ddply(xrib11,c("Year","SSP"),transform,v80=quantile(Value,0.8))
      xrib13 <- ddply(xrib12,c("Year","SSP"),transform,v90=quantile(Value,0.9))
      xrib14 <- ddply(xrib13,c("Year","SSP"),transform,v95=quantile(Value,0.95))
      xrib15 <- ddply(xrib14,c("Year","SSP"),transform,v98=quantile(Value,0.98))
      xrib16 <- ddply(xrib15,c("Year","SSP"),transform,v99=quantile(Value,0.99))
      
      names(xrib16)[8]<-paste("max")
      names(xrib16)[9]<-paste("min")
      
      #            col<- c('#00FF00', '#0000FF',"#000000")
      col<- c('mediumblue','orangered3',"#000000")

      xrib16[1] <- lapply(xrib16[1],as.integer)
      xrib16[1] <- (xrib16[1]-18)*10+2010
      xrib17 <- xrib16 %>% 
        distinct(Year,SSP,GCM,XXX,max,min,v1,v2,v5,v10,v20,v33,v50,v66,v80,v90,v95,v98,v99)
      #2015 continuation
      fao2015 <- filter(faoworld,Year==2010)
      fao2015 <- as.numeric(fao2015[1,3])
      xrib16add <- xrib17 %>%
        mutate(v000=fao2015,v00=fao2015,v01=fao2015,v1=fao2015,v0=fao2015,v2=fao2015,v3=fao2015,v4=fao2015,v5=fao2015,v6=fao2015,v7=fao2015,v8=fao2015,v9=fao2015,v10=fao2015,v11=fao2015,v12=fao2015,v13=fao2015,v14=fao2015,v15=fao2015) %>%
        select(-max,-min,-v1,-v2,-v5,-v10,-v20,-v33,-v50,-v66,-v80,-v90,-v95,-v98,-v99)
      names(xrib16add) <- c("Year","SSP","GCM","XXX","max","min","v1","v2","v5","v10","v20","v33","v50","v66","v80","v90","v95","v98","v99")
      xrib16add[1] <- 2010
      xrib18 <- rbind(xrib17,xrib16add)
      plotrib <-ggplot()+ 
#        geom_ribbon(data=xrib18,aes(x=Year,ymin=min, ymax=max, fill = SSP), alpha=0.1) +
        geom_ribbon(data=xrib18,aes(x=Year,ymin=v1, ymax=v99, fill = SSP), alpha=0.2) +
        geom_ribbon(data=xrib18,aes(x=Year,ymin=v5, ymax=v95, fill = SSP), alpha=0.3) +
        geom_ribbon(data=xrib18,aes(x=Year,ymin=v10, ymax=v90, fill = SSP), alpha=0.4) +
        geom_ribbon(data=xrib18,aes(x=Year,ymin=v33, ymax=v66, fill = SSP), alpha=0.7) +
        geom_line(data=xrib18,aes(x=Year,y=v50,group = SSP, colour=SSP), size = 2) +
        labs(aes(colour = SSP)) + 
        scale_colour_manual(name="Median value",values=col) + 
        scale_fill_manual(name="Impact range",values=col) +
        theme_bw(base_size = 30) + labs(y =labelname, x="") +
#              ylim(minv-abs(minv)*0.01, maxv+abs(maxv)*0.01) +
        MyThemeLine +
#        ggtitle(mainlabel) +
        geom_line(data=faoworld2,aes(x=Year,y=Value,group=Reg,xend=2010)) 
#        xlim(1990,2050)
#        ylim(minv-abs(minv)*0.01, maxv+abs(maxv)*0.01) 
      
      plot(plotrib)
      
      filenamehist <- paste('../output/png/', varname ,'/',regionname,'_',yearname,cropname,'ribbon_rcp85.png', sep = '')
      ggsave(plotrib,filename=filenamehist, dpi=250,width=8, height=6)
      #        }
    }
  }
}
#----return period
var.in <- c("PCAL")
varname <- var.in
y <- pdata[pdata$VAR==varname & pdata$SSP=="SSP2" & pdata$Crop=="dum", c(-1,-4,-5,-8,-9,-10)]

names(y) <- c("Year","Region","RCP","CO2fert","Value")	#Rename the column
interact<-interaction(y$RCP, y$CO2fert, sep=" : ")
yall <- 0

returnp<-c(0.01,0.02,0.05,0.1,0.2,0.33)
returnn<-c("RP001","RP002","RP005","RP01","RP02","RP033")
for(i in 1 : 6){
  a <- returnp[i]
  b <- returnn[i]
  y5 <- ddply(y,c("Region","Year","interact"),transform,b,v1=quantile(Value,a))
  names(y5) <- c("Interact","Year","Region","RCP","CO2fert","Value","V","return")	#Rename the column
  if (yall==0){
    yall <- y5
  }else{
    yall <- rbind(yall,y5) 
  }
  y5<-0
}
zall <- yall[c(-1,-6)]
z <- unique(zall)
symDim <- 6
attr(z, "symName") <- "Foodshock"
lst <- wgdx.reshape(z, symDim)
wgdx.reshape(z, symDim, gdxName = "../data/FoodshockR.gdx")
system(paste("gams","../prog/foodstock.gms",sep=" "))

linepalette <- c("#4DAF4A","#377EB8","#E41A1C","#FF7F00","#984EA3","#FFFF33","#A65628","#F781BF","#8DD3C7","#FB8072","#80B1D3","#FDB462","#B3DE69","#FCCDE5","#D9D9D9","#BC80BD","#CCEBC5","#FFED6F")
ReturnPeriodFS <- rgdx.param('../output/gdx/foodstock.gdx','foodstockEJ')
ReturnPeriodpr <- rgdx.param('../output/gdx/foodstock.gdx','foodstockR_per')
ReturnPeriodpr <- ReturnPeriodpr[c(7)]
#ReturnPeriodFS[6] <- ReturnPeriodFS[6]*100
ReturnPeriod <- cbind(ReturnPeriodFS,ReturnPeriodpr)

x <- ReturnPeriod[ReturnPeriod$Y=="2050" & ReturnPeriod$co2fert=="noco2" & ReturnPeriod$SSP=="SSP3" ,c(-1,-4,-5,-6)]
names(x) <- c("Region","RCP","Stock","RetPer")
x2 <- x[x$Region=="WLD",c(2,3,4)]
plot1 <- ggplot() + 
  geom_point(data=x, aes(x=RetPer,y=Stock,group=RCP,color=RCP)) + MyThemeLine +
  xlab("return period (Year)")+ylab("Required food stock (proportion to world wheat production(%)")
plot2 <- plot1+ facet_wrap( ~ Region,scales="free")

ggsave(plot2,filename="../output/png/return_stock.png", dpi=250,width=10, height=6)

plot3 <- ggplot() + 
  geom_bar(data=x[x$Region!="WLD" & x$RCP=="RCP85",c(1,2,3,4)], aes(x=RetPer,y=Stock,fill = Region),stat = "identity") + MyThemeLine + scale_fill_manual(values =linepalette)+
  xlab("return period (Year)")+
#  ylab("Required food stock (proportion to world wheat production(%)") +
  ylab("Required food stock (EJ)") +
  geom_line(data=x[x$Region=="WLD" & x$RCP=="RCP85",c(1,2,3,4)], aes(x=RetPer,y=Stock,color=Region),stat = "identity") + MyThemeLine + scale_color_manual(values =c("red")) +
  annotate("segment",x=0,xend=100,y=0,yend=0,linetype="solid",color="grey")

ggsave(plot3,filename="../output/png/return_stock_reg.png", dpi=250,width=10, height=6)

#---end return period


#----- PPOP_UNSH & PCAL: Boxplot across regions ----
var.in <- c("PCAL")
var.in <- c("PPOP_UNSH")
var.in <- c("PCAL","PPOP_UNSH")

for(k in 1:length(var.in)){
  varname <- var.in[k]
  
  
  if (varname=="PCAL"){
    labelname <- c("calorie intake [kcal/person/day]")
  } else if (varname=="PPOP_UNSH"){
    labelname <- c("population at risk of hunger [billion]")
  } else if (varname=="YIELD"){
    labelname <- c("yield [tonne/ha]")
  } else {
    labelname <- c("aaaaaaaa")
  }
  
  filenamebox <- paste('../output/png/', varname ,'/',yearname,'.png', sep = '')
  y <- pdata[pdata$VAR==varname & pdata$Year==yearname  & pdata$SSP=="SSP2" & pdata$Crop=="dum" & pdata$Region!="WLD", c(-1,-2,-4,-10)]
  
  names(y) <- c("Region","GCM","RCP","CO2fert","SYR","N","Value","NoCC")	#Rename the column
  interact<-interaction(y$RCP, y$CO2fert, sep=" : ")
  
  minv <- min(y[7],y[8])
  maxv <- max(y[7],y[8])
  
  #	bw <- diff(range(y$Value))/50
  #	mu <- ddply(y, "RCP", summarise, grp.mean=mean(Value))
  #	mu1 <- ddply(ypre, "RCP", summarise, grp1.mean=mean(Value))
  
  
  y0 <- ddply(y,c("Region","interact"),transform,v0=min(Value))
  y1 <- ddply(y0,c("Region","interact"),transform,v25=quantile(Value,0.25))
  y2 <- ddply(y1,c("Region","interact"),transform,v50=median(Value))
  y3 <- ddply(y2,c("Region","interact"),transform,v75=quantile(Value,0.75))
  y4 <- ddply(y3,c("Region","interact"),transform,v100=max(Value))
  z <- regionlabel(y4)
  
  gbox <- ggplot() + 
    geom_boxplot(data=z, aes(x=Region,y=Value, color=interact, ymin = v0, lower = v25, middle = v50, upper = v75, ymax = v100), fill="grey80", stat="identity") +
    ylim(minv-abs(minv)*0.01, maxv+abs(maxv)*0.01) +
    #	ggtitle(mainlabel)
    ylab(labelname) + xlab("") +	MyThemeLine +
    geom_point(data=z, aes(x=Region, y=NoCC),size=10,shape=73,colour="black")
  
  gbox1 <- gbox + coord_flip()
  #  g2 <- g1 + facet_wrap(~RCP,ncol=2)
  
  plot(gbox1)
  
  ggsave(gbox1,filename=filenamebox, dpi=250,width=10, height=6)
  
}

#----End --- PPOP_UNSH & PCAL: Boxplot across regions ----

#-----  YIELD(YEXO) & price(XPRP): Boxplot across regions ----


var.in <- c("XPRP")
var.in <- c("YEXO")

for(k in 1:length(var.in)){
  varname <- var.in[k]
  
  if (varname=="PCAL"){
    labelname <- c("calorie intake [kcal/person/day]")
  } else if (varname=="PPOP_UNSH"){
    labelname <- c("population at risk of hunger [billion]")
  } else if (varname=="YEXO"){
    labelname <- c("yield [tonne/ha]")
  } else if (varname=="XPRP"){
    labelname <- c("production price [2005=1]")
  } else {
    labelname <- c("aaaaaaaa")
  }
  
  #  filenamebox2 <- paste('../output/png/', varname ,'/',yearname,'_limited.png', sep = '')
  filenamebox2 <- paste('../output/png/', varname ,'/',yearname,'.png', sep = '')
  
  #colnames <- c("VAR","Year","Region","SSP","GCM","RCP","CO2fert","SYR","N","Crop","Value","NoCC")
  
  # Limit the num of regions
  #pyld <- pdata[pdata$VAR=="YEXO" & pdata$Year==yearname  & pdata$SSP=="SSP2" & pdata$Region!="WLD" & 
  #               pdata$Region!="XOC" & pdata$Region!="XNF" & pdata$Region!="XER" & pdata$Region!="XE25" & 
  #               pdata$Region!="USA" & pdata$Region!="TUR" & pdata$Region!="JPN" & pdata$Region!="CIS" & 
  #               pdata$Region!="CAN"& pdata$Region!="BRA" & pdata$Region!="CHN", c(-1,-2,-4)]
  
  # Limit the num of regions
  pyld <- pdata[pdata$VAR==varname & pdata$Year==yearname  & pdata$SSP=="SSP2" & pdata$Region!="WLD", c(-1,-2,-4)]
  
  
  names(pyld) <- c("Region","GCM","RCP","CO2fert","SYR","N","Crop","Value","NoCC")	#Rename the column
  
  minv <- min(pyld[8],pyld[9])
  maxv <- max(pyld[8],pyld[9])
  #  maxv <- c(20)
  
  #	bw <- diff(range(y$Value))/50
  #	mu <- ddply(y, "RCP", summarise, grp.mean=mean(Value))
  #	mu1 <- ddply(ypre, "RCP", summarise, grp1.mean=mean(Value))
  
  interact<-interaction(pyld$RCP, pyld$CO2fert, sep=" : ")
  
  pyld0 <- ddply(pyld,c("Region","interact","Crop"),transform,v0=min(Value))
  pyld1 <- ddply(pyld0,c("Region","interact","Crop"),transform,v25=quantile(Value,0.25))
  pyld2 <- ddply(pyld1,c("Region","interact","Crop"),transform,v50=median(Value))
  pyld3 <- ddply(pyld2,c("Region","interact","Crop"),transform,v75=quantile(Value,0.75))
  pyld4 <- ddply(pyld3,c("Region","interact","Crop"),transform,v100=max(Value))
  pyld5 <- regionlabel(pyld4)
  pyldz <- croplabel(pyld5)
  
  
  gyldbox <- ggplot() + 
    geom_boxplot(data=pyldz, aes(x=Crop, y=Value, color=interact, ymin = v0, lower = v25, middle = v50, upper = v75, ymax = v100), fill="grey80", stat="identity") +
    ylim(minv-abs(minv)*0.01, maxv+abs(maxv)*0.01) + ggtitle("") +
    ylab(labelname) + xlab("") +	MyThemeLine +
    geom_point(data=pyldz, aes(x=Crop, y=NoCC),size=6,shape=73,colour="black")
  
  gyldbox1 <- gyldbox + coord_flip()
  gyldbox2 <- gyldbox1 + facet_wrap(~Region, ncol=4)
  
  plot(gyldbox2)
  
  ggsave(gyldbox2,filename=filenamebox2, dpi=250,width=7.5, height=7.5)
  
}



#----- END --- YIELD : Boxplot across regions ----






